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SUMMARY 

2 i 

A method is presented for the calculation of multicomponent plasma 
properties. The method includes both Debye-Buckel and ionization potential 
lovering effects . These effects are shovn to be partially compensating. A 
perturbation technique is used to make the calculation compatible vith a 
scheme used for calculating ideal gas properties. Thermodynamic derivatives, 
such as heat capacity, can be calculated directly. 

INTRODUCTION 

Plasmas have become important in many areas of scientific and technical 
vork. Often the important plasma systems are not simple one-component systems 
but may involve many different chemical elements. There axe several schemes 
(ref. l) and many computer programs for the computation of the equilibrium 
properties of thermodynamic systems that can be regarded as mixtures of ideal 
gases. These methods, however, must be modified to account for the presence of 
Coulomb forces. 

In principle, the presence of long-range Coulomb forces in a plasma in- 
validates the independent-particle (ideal gas) approach and requires the con- 
sideration of interactions among all particles in the system. Even for simple 
systems, hovever, this approach is extremely difficult and has yielded an 
exact solution only for a one-dimensional plasma model (refs. 2 and 3). A 
scheme is presented in this report for calculating the thermodynamic properties 
of many component plasmas that can be applied to most systems and can be ex- 
pected to yield results vhich are an improvement over the ideal gas calcula- 
tions . 

The thermodynamic properties for a system of independent particles (ideal 
gas) can be calculated from a paxtition function that is a product of tvo 
factors. One factor represents the contribution of the translational 

kinetic energy to the partition function vhile the other factor Qjn-h repre- 
sents the contributions from the internal energy levels . Because of the 
assumed independence of the particles, the energy levels used to calculate the 
internal partition functions are the energy levels for an isolated particle. 
These energy levels are available for many species from spectroscopic data. 
Plasma computations that assume ideal gas behavior and neglect Coulomb effects 





are typified by calculations of Kubin and Presley (ref. 4) and 'early computa- 
tions of Hilsenrath (ref. 5) and Rouse (refs. 6 and 7). A calculation that 
assumes nonideal behavior but still neglects Coulomb effects is reported by 
Rosenbaum and Levitt (ref. 8). 

An attempt to incorporate the effect of Coulomb forces produces two changes 
in the partition function from the partition function in the ideal gas case. 
First, a third factor appears in the partition function. This factor accounts 
for the effect of interparticle forces upon the translational partition function 
and represents the contribution to the partition function of the classical con- 
figuration integral. Second, in calculating the internal partition function, 
one should no longer use the energy levels of an isolated atom or molecule but 
should use the energy levels of the atom or molecule in the external field pro- 
duced by the other particles of the system (refs. 9, 10, and ll). Thus, for a 
plasma, the internal partition function is dependent upon the amount of ioniza- 
tion in the plasma and therefore is a function of the charged particle concen- 
trations. The concentration dependence of the energy levels is generally un- 
known. Some recent attempts have been made to estimate this effect for 
hydrogen-like atoms by calculating the energy levels for a screened Coulomb 
potential (refs. 10 and 12 to 15). No data of this type, however, exist for 
other species. An alternate and somewhat simpler technique has also been used 
to take into account the fact that the internal partition function has concen- 
tration dependence. This method uses isolated particle energy levels together 
with an energy level cutoff and an ionization potential lowering that are 
dependent upon electron concentration; it has been used by Rouse (refs. 16 
to 18), Drellishak, Knopp, and Cambel (ref. 19), and Drellishak (ref. 20). The 
difference between the two procedures, then, is that the first adjusts internal 
energy levels by varying amounts while the second, in effect, adjust the levels 
by a constant amount. When only the ground state is used, as by Harris 
(ref. 21 ), the two procedures are equivalent. 

The previous discussion was concerned with the effect of Coulomb forces on 
the internal partition function Q -j n f. . The effect of these forces on the 
translational partition function Q^ r is usually calculated by assuming the 
validity of the Debye-Huckel approximation. The Debye-Huckel theory is dis- 
cussed in various textbooks (e.g., refs. 22 and 23). Its validity is discussed 
by Balazs (ref. 24) in connection with the one -dimensional plasma for which an 
exact solution is available. Duclos and Cambel (ref. 25) summarize the con- 
clusions of various authors regarding the region of validity of the Debye-Huckel 
approximation for three-dimensional plasmas. Myers, Buss, and Benson (ref. 26), 
Hilsenrath and Klein (ref. 27), Harris (ref. 28), Harris and Trulio (ref. 29), 
and McGee and Heller (refs. 30 and 3l), among others, have all used the Debye- 
Huckel approximation in calculations . Unfortunately, the calculations of McGee 
and Heller are in error since they mistakenly use the excess Helmholtz free 
energy for the excess Gibbs free energy and, in addition, they use an incorrect 
expression for the excess entropy. 

Some confusion exists in the literature regarding the term "ionization 
potential lowering." This term was discussed by Ecker and Kroll (ref. 32), 
who showed that both the translational and internal energy level effects could 
be regarded as an ionization potential lowering. This was also mentioned by 
Harris (ref. 2l). In this report the term "ionization potential lowering" will 


2 



mean only the effects associated with internal energy levels . 

! The previously discussed Coulomb effects on both and Qint are in- 

cluded in this report. For the effect on fya*, validity of the Debye-Huckel 
limiting law approximation is assumed. For the effect on Qj. n t> an ionization 
potential lowering that is proportional to k , the reciprocal of the Debye 
! length, is assumed. This assumption is equivalent to that of Harris (ref. 2l) 
and Drellishak, Khopp, and Cambel (ref. 19). In performing numerical computa- 
tions, some scheme must also be used for cutting off the infinite number of 
unperturbed internal energy levels. This aspect of the problem will not be 
considered further; however, several of these cutoff schemes are reviewed by 
j Hester and Sewell (ref. 33) and also by Margenau and Lewis (ref. 11 ). 

The usual procedure for calculating the thermodynamic properties of plasmas 
is to assign volume and temperature as the independent thermodynamic parameters. 
In contrast, the method developed in this report assigns the pressure and one 
, other thermodynamic variable selected from the three quantities temperature, 
enthalpy, and entropy as independent parameters. For many applications, it is 
more convenient to have data as functions of pressure rather than volume. In 
addition, this procedure permits the calculation of thermodynamic derivatives 
such as specific heat without resorting to numerical differentiations of tabu- 
: lar data (ref. 19). 

THERMODYNAMIC JUNCTIONS FOR MULTICOMPONENT PLASMA 

The Helmholtz free energy A of a plasma may be written as 

! A = A-j- + A c (l) 

where the subscript I denotes the ideal gas contribution and c denotes the 
contribution from Coulomb forces. (All symbols are defined in the appendix.) 

In terms of the canonical partition function for an ideal gas system, 

Aj. = (2) 


where T is the absolute temperature and k is Boltzmann's constant. The 
translational partition function is given by 



(3) 


where V is the volume, ft is Planck's constant divided by 2 jt, m^ is the 
mass per particle of species i, and % is the number of particles of 
species i. The internal partition function is given by 


Qint ~ 



where the internal partition function for the 


^th 


species is 


U) 


3 



(5) 
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where Eq^ ^ is the ground state energy for the i th isolated species and 

are the isolated energy levels relative to the ground state; Qj.nt is inde- 
pendent of volume. 


The contribution of the Coulomb forces to the Helmholtz free energy can be 
written in the form 



( 6 ) 


i 

The first term in this equation represents the contribution of the Debye-Huckel 
limiting law (refs. 22 and 23); that is, it is the effect of Coulomb forces on 
the translational partition function. The second term represents the contribu- 
tion from the displacement of the internal energy levels. It arises because the 
truncated partition function goes to exp( -AE^ /kT) when each in- 

ternal energy level is displaced by the same amount AE (•*•). The amount of 
energy level displacement will be a function of ionization. The quantity K is 
the reciprocal of the Debye length and is given by 


K 


? 4:rt V 
” V>T 


VkT 


N ± q£ 


(7) 


where qj_ is the charge on the i th species. The dielectric constant that 
usually appears in the definition of the Debye length has been set equal to 
unity in equation (7). The displacement of the internal energy levels of the 
i^k species is assumed to be proportional to k, and the proportionality con- 
stants Pj_ may be different for each species: 

AE (i ) = p iK a 0 (8) 


where has the dimensions of energy and a^ is the Bohr radius . Combining 

equations (6) and (8) gives 





Pi W i 


( 9 ) 


The assumption (eq. (8)) is consistent with the assumptions made by 
Harris (ref. 2l) and Drellishak, Knopp, and Cambel (ref. 19). For the hydrogen 
plasma, Harris used |3 jj = 2Ijj where Ijj is the hydrogen atom ionization po- 
tential. For species in the argon plasma, Drellishak assumed that the change 
in ionization potential AI was given by the expression 


4 


3 — 0,lj2,3, 4 


( 10 ) 







If it is assumed that the energy of the free electron is zero, the ionization 
potential of successively ionized species of argon can he expressed in terms 
of their isolated particle ground state energies 

!(Ar+j) = “ 4 J) <5 = °A> 2 > 3 > 4 

Therefore, the change in ionization potential is 


Ai = ae(J +1 ) - ae(^ 

(Ar + J ) 


j = 0,1,2, 3, 4 


The use of equations (8), (10), and (ll) implies the following p's: 
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( 11 ) 


( 12 ) 


P (Ar +5 ) ° 

Once the form of the Helmholtz free energy (eg.. (9)) has been selected, all 
other thermodynamic functions can be calculated from it by the usual thermo- 
dynamic relations. For example, with U = -T 2 [S(A/T)/ST] V the internal 
energy can be calculated as * 1 

U = Up + U c (13) 


where 


kT 


8 * 


Vk" 


+ 3 
2 


Ka 0 

kT 


L Pi N i 


(14) 


The equation of state for the system is obtained from the thermodynamic 
relation P = -(dA./dV)j^if and can be written in the form 

I>i N i 

PV _ „ „ 1 Vk 3 . 1 Ka 0 i 


W?E Ni 

i 


s Z = 1 - 


24*£il + 2 kT ^ Ni 


(15) 


The deviation of Z from unity indicates the departure from ideal behavior due 
to Coulomb forces, and thus Z will be called the Coulomb compressibility. 

The two terms in equation (15) that account for the deviation from unity are 
compensat ing if the p^ axe taken to be positive. In terms of Z, equa- 
tions (9) and (14) become 
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( 16 ) 


4 



1) XX 

i 





i) 2X 

i 


(17) 


The enthalpy H and the ideal gas enthalpy Hj are defined as 

H=U+FV=U I +U C +PV 

Hj = Uj + kT^X 

i 


Equations (15), (17), (18), and (19) may he combined to give 

H = Hj + 4(Z - l)kTXX 

i 


(18) 

(19) 


( 20 ) 


The expression for entropy in terms of volume and temperature can be ob- 
tained from the thermodynamic relation S = -(dA/dT)y jj.. Just as in the case 

of the other thermodynamic functions, the contribution of the Coulomb forces can 
be written in terms of the Coulomb compressibility: 

S = Sj + (Z - 1 + In Z)k £X (21) 

i 

The In Z term in the correction for nonideality arises because the volume was 
eliminated by means of the equations of state (15). Usually In Z is ignored; 
however, by doing this one is disregarding a term that is of the same order of 
magnitude as the term retained. 

Expressing the equations for enthalpy and entropy in terms of moles n^ 
rather than the number of particles Np gives 

m m 

H = Yj M^i + 4 ( Z - 1 ) RT Yj ni (22) 

i=l i=l 

m m 

s = Y ( ®T ) i n i + (z - 1 + in Z)K Y n i (23) 

i=l i=l 

where 


(St) ± = - R In pi 


6 


n i p 

p, = _sL_ 
j in 


n-; 


i=l 


In these equations it has been assumed that there are m different species 
in the system. With the use of equations (22) and (23) the Gibbs free energy 
becomes 

F = Fp + F c (24) 

where 

p i = t 

i=l 

m 

F c = [3(Z - 1) - In z]RT X) n i (26) 

i=l 


INVERSION OF EQUATION OF STATE 

It is desired to express the equations of this report in terms of pres- 
sure, rather than in terms of volume, as one of the independent parameters. 

Thus an expression for volume in terms of pressure is necessary in order to 
eliminate volume from all thermodynamic functions. This expression can be ob- 
tained by inverting the equation of state. This inversion is equivalent to ex- 
pressing the Coulomb compressibility in terms of pressure and temperature. 

If the dimensionless ionization parameter 



is introduced, the equation of state (15 ) becomes 

1 + oZ -1 / 2 = Z (28) 

The analytical solution of equation (28) in the form Z = Z(a) can be found 
conveniently by studying the associated cubic equation 

Z 3 - 2Z 2 + Z - a 2 * 0 (29) 

This equation was obtained from equation (28) by rearranging terms and then 
squaring. Not all roots of the cubic equation (29) are valid solutions for the 
equation of state ( 28 ) . The solution that is valid for both equations is 


7 




where 

2 

cos cp = -1 + 

a 0 

where 

a ° = 

Figure 1 is a plot of Z as a function of a. If the ionization poten- 
tial lowering is neglected, a < 0 and Z < 1; whereas, if the Debye -Huckel 
limiting law is neglected, a > 0 and Z > 1. Because these two contributions 
are partially compensating, considerable ionization can occur with little effect 
on the equation of state. It may be noted that there is no solution for 
a < -ao* At a = -a q, Z has the value l/3 while, at a = cgq> Z equals 4/3. 

For | a | < 0.1, the trigonometric solution (eq. (30)) can be replaced by 
the expansion 


Z = 1 + a 


+ 1 3 _ 
2 8 


a 


231 a 5 
128 


I a, I £ 0.1 


This expansion gives at least seven figure accuracy in this interval. 


(31) 


ITERATION EQUATIONS 

In this section, the iteration equations used to calculate chemical equi- 
librium will be derived. In the following section, the equations required for 
the calculation of thermodynamic derivatives will be given. For the purposes 
of this report, the derivations will be based on the minimization of the Gibbs 
free energy and will be related to the somewhat different free energy minimiza- 
tion methods described by White, Johnson, and Dantzig (ref. 34) and Zeleznik 
and Gordon (ref. l). In order to make the method compatible with the tech- 
niques used in the computer programs described by Gordon and Zeleznik 
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Coulomb compressibility, Z 



(refs. 35 and 36), the derivation used here 
will be such that in the ideal gas limit the 
equations will be identical in form to those 
obtained by the equilibrium constant method 
(ref. l). This can be accomplished by re- 
garding the sum of the moles of all gaseous 
species as fixed at a value that is numeri- 
cally equal to the assigned pressure Pq. 

This assumption requires treating the total 
mass of reactants as a variable just as has 
been done in the past for the ideal gas 
(refs. 1 and 37). Since the derivation of 
the iteration equations closely parallels the 
minimization procedure of reference 1, many 
of the details will be omitted. 


Figure L - Coulomb compressibility as function 

of ionization parameter. The thermodynamic state is assumed to be 

specified by assigning the pressure and one 
other thermody nami c variable selected from the three quantities, temperature, 
enthalpy, and entropy. The pressure of the system is specified by 


P Q - P = /JsP = 0 


(32) 


m 

where P = 53 Pi* ®h e specification of the thermodynamic state is completed 
i=l 

by assigning either temperature, enthalpy, or entropy as follows: 


T 0 - T = 0 
hQ - h = = 0 

Sq — s = - 0 

where 

m 

Ah. £ 

i=l 

m m 

As = (St)-^ + (Z - 1 + In Z)R ^ n^ 

i=l fcl 



(33) 

(34) 

(35) 

(36) 

(37) 


The quantity A is the total mass of the system. Since the size of the system 
is arbitrary, it is convenient to let A be that mass which is required in 
order that the partial pressure of each species is numerically equal to the 
number of moles of that species ; that is, pi = ni-_ In the following equa- 
tions, nj_ is substituted for pj_ . The value of A that permits this sub- 
stitution is determined simultaneously with the other variables of the problem. 


Once the thermodynamic state has been specified, the equilibrium composi- 
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tions can be obtained by minimizing the Gibbs free energy subject, to mass - 
balance constraints. If there are (Z - l) different chemical elements (where 
Z is some fixed integer), the Z equations for conservation of species and 
charge neutrality can be written as 

bj - b ± = &> ± = 0 (i = 1 . . . Z) (38) 

m 

where Abj_ = ^ a^nj . The a^j are the stoichiometric coefficients. In 

the set of equations (38), the equation corresponding to i = Z refers to the 
conservation of charge. For the case of a neutral plasma, bf = 0. 

If the state is specified by assigning a temperature and a pressure, the 
iteration equations are obtained by minimizing the quadratic approximation to 
F/RT subject to the constraints (eqs. (32) and (33)) and the linearized form 
of equation (38). The mass balance constraint is incorporated by means of 
Lagrange multipliers whereas the pressure equation and the temperature 

equation are used directly. The resulting iteration equations are 


m Z 

E + r kj) A in n j - E *i a ik 
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It may be noted that, for assigned temperature, equation (42) is trivial and 
may be omitted together with the A In T term in equation (39). If, however, 
the temperature is specified indirectly by means of equation (34), then the 
linearized form of equation (34) is used in place of equation (42): 



The following equation, (45), is used when the temperature is specified in- 
directly by means of equation (35) instead of equation (42). This equation is 
the linearized form of equation (35) to which has been added equation (4l) 
multiplied by the factor (2 - Z - In Z): 


Z 


r 

n j 

M 

+ nj(2 - Z - In Z) ► 

T,P 



J 


s 

" r 


+ T 


Vs 

\R 


"St 




A In T = —5— + AP( 2 - Z - In Z) (45) 


The most important difference between the preceding sets of iteration 
equations and the corresponding equations for the case of an ideal gas is the 
appearance of the matrix T in equation (39). In the ideal case (P = 0), 
these equations could be used to eliminate A In n^ ( j = 1 . . . m) in terms of 
(i = 1 . . . 7), A In A, and A In T. Thus the number of equations could be 
reduced from m+Z + 2 to 1 + 2 and the amount of numerical computation 
could be considerably reduced. For the case of a plasma (r ^ 0), such a 
reduction cannot be performed as easily. This problem will be discussed more 
fully in a later section. 


THERMODmMIC FIRST DERIVATIVES 

In addition to the thermodynamic functions themselves, it is often 
desirable to have available the first derivatives of the thermodynamic func- 
tions. All possible first derivatives can be expressed in terms of any three 
independent ones. As in earlier papers (refs. 35 and 37), the heat capacity at 
constant pressure 
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= /° H \ _ /dH\ 
P \5T/p " \^T/p 


SH \ 
,^T,P 


a In n^ 


and the derivatives of volume with respect to temperature and pressure ( cN /c)T )p 
and (Sv/^P)ip are selected as the independent derivatives. If an average 
molecular weight is introduced, the equations of state can be written 


PV = Z^RT 


This equation leads to the relations 


a In V \ „ 1 _ ( a In M \ + 1 dZ 


In T/p Z da 


(^)p,n 1 


ILL 

z- 


a In n ± 


a In 
,a In Pi 


\ = . r + / d in m\ 1 + I az 1 7 act \ 
y T I y a In P /<pj Z da |\a In Py 


V n (to \ / a »i\ ' 

/ “Assy u to p i 

t X X/ T,r 7 tJ 


Prom the procedure of reference (37), it is readily established that 

/a In m\ /a In a\ 


o In n 


d In A 





Therefore, the evaluation of the three independent first derivatives reduces 
to the problem of calculating derivatives with respect to In T at constant 
pressure -and with respect to In A at constant temperature. The derivatives 
with respect to In T are obtained by solving the set of equations 

m l 

Z au k / a m nA au k ( a» ± \ ( au k \ _ 

n j drij yd In T /p ? \d In T/p yd In T/p 


3=1 


m 




(k = 1,2,3 . . . m) 
(i = 1,2,3 . . . l) 


3=1 


m 




= 0 


3=1 


where 


Uk ■ db 


fl + E n / Ab i -^ a iJ D 3 

\ 3=1 



Vet/ 

) 

£tak 


E 




T,P i=l 


T,P 


(53) 


(54) 


(55) 


(56) 


These m + l + 1 equations are identical in form to the m + l + 1 iteration 
equations (39) to (41) although they are not derived from the iteration equa- 
tions. They can be thought of, however, as being obtained from the iteration 
equations by setting the right-hand side equal to zero and by formally dividing 
through by A In T. 


The derivatives with respect to In A are obtained by solving the set of 
equations 




(k = 1,2,3 . . 


. m) 


(57) 



3=1 


(i = 1,2,3 . . . 1) 


(58) 
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Again there is an obvious similarity between these equations and the iteration 
equations ( 39 ) and ( 40 ) . 


SOLUTION OF EQUATIONS BY PERTURBATION 

The problem of calculating the thermodynamic properties of a multicomponent 
plasma has been reduced to the problem of solving one of three different sets 
of simultaneous linear iteration equations. The set of equations (39), (40), 

( 4l), (42) is used if the thermodynamic state is specified by assigning a 
temperature and a pressure. A second set (eqs. (39), (40), (4l), (44)) is used 
if the state is specified by enthalpy and pressure while the third set (eqs. 
(39), (40), (4l), (45)) is used when the state is specified by entropy and 
pressure. When the iteration procedure has converged, two additional sets of 
linear equations must be solved in order to calculate the three independent 
derivatives discussed earlier. One of these sets is composed of equations (53), 
(54), (55) while the other is formed from equations (57), (58). 

These five sets of linear equations are very similar to each other, and 
therefore simultaneous discussion of the solution of these sets is possible. 

In the discussion of these equations, the more concise matrix notation 
M v = w will be used to symbolize any one of the five sets of linear equations . 
Further the matrix M and the column vectors v and w will be assumed to 
have been partitioned so that the set of equations can be written in the form 


fh.1 _ /-i\ 

V^ivW W 


The col umn vector v^ is associated with the variables related to 
n<(j = 1,2, . . .,m), while v is associated with the remaining variables. 

The matrix 7^ is 

square and has the dimensions 7 + 2 for the iteration equations, while for 
the sets (53), (54), (55) and (57), (58) it has the dimensions 7+1 and 7, 
respectively. 

If the set of equations represented by equation (59) is relatively small, 
one can work directly with these equations. In multicomponent systems, however, 
the number of equations could be quite large and their solutions could involve 
a considerable amount of numerical difficulty. For this reason, working with 
a reduced set of equations analogous to those used for the ideal gas problem 
would be convenient. As has already been pointed out, the principal difference 
between equation (59) and the corresponding ideal gas equations lies in the 
fact in the ideal gas case was a unit matrix. This permitted the first 

member of equation (59) to be used in the direct elimination of v^ from the 
second member of the set and thereby gave a considerably reduced set of linear 
equations in the variables v ^ • Since the objective is to provide a calculating 
scheme compatible with the ideal gas case, a reduced set of linear equations 


square and has the dimension m. The matrix VL ? is also 
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must first be 'constructed. With the square submatrix JL assumed to be non- 
singular ‘ LL 



(60) 


This equation is used to eliminate v from the second member of equation (59), 
which yields ^ 

(—22 " ^ I ^ LI ^ Le ) -2 = -2 ^ 1 ^ 1-1 

Equations (6l) possess the disadvantage that they contain the matrix 
= (I + r) -1 . The need for numerically calculating M£^- would be 

eliminated if a closed-form expression could be obtained for M^. Obtaining 
this expression is certainly not possible in a general case, bu,t it can be done 
when the contribution of T to M can be regarded as a perturbation . Then 

the iterative expansion of M -x is 



(62) 


where r° = I. 


If the notation N v = is used for the reduced set of equations (6l) 

C* -» 

and if equation (62) is vised to eliminate M~^ from both N and y, 


W = ^ 


-\ 


k=0 


= £ 


k=0 


where 


= M - ( -l) k M 

- -22 -21- “L2 

- (-i) k *L,r k w, 

*- -2 -21 1 




J 


(63) 


(64) 


The index k effectively gives the order to which the perturbation matrix T 
appears in the various terms. Therefore, conventional perturbation theory can 
now be applied in the solution of the reduced set of equations (6l). Writing 
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v = y v (k) 

~ 2 U - 2 


( 65 ) 


and equating terms in equations (6l) with the same order of perturbation yields 
the usual equations 


n(0)v(°) = Z (°) 


k = 0 


M (0) y £ , 


s k> - z (k> - Z S U) l2 k '' i) 


k = 1,2,3 . • . 


(66) 




These equations can be used to calculate to the desired order of perturba- 
tion. An expression for can then be obtained from by using equa- 

tion (60). Grouping terms with the same order gives 


00 , N 00 
n=0 n=0 L 


(-1 ) n r n w 1 - ^ (-i) k ^ 2 v^ n - k ) 


k=0 


(67) 


The explicit expressions for the matrices N^ k ) and the vectors k will 
be given for the case in which the equations to be solved by perturbation 
correspond to the iteration equations . Using the notation 


a k = 


\ da ’ 

K 


= a I 


2 

^k 


T,P 


m 


Zniq? 


3 

2 m 


q k - 3a 0 p k 


i=l 


mm m 

Y n i Y n i<il - 3a 0 y p ini 

i=l i=l i=l y 


( 68 ) 


5a 




= -2a 


(69) 


yjLY 

°\RT/ 


IS 


l k 


Jt,p 


( F t ) 1 


RT 


- * - 1 


+ hin, + 1 + 3(Z - 1) - hi Z 


dZ 


m 

+ * Z Yj ni da ak 

i=l 


(5Z - 1) 


(70) 


(71) 
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h k = 


<£) 


dn-b 


m 


JT,P 


RT 


+ 4(Z-l) + 4£n i |a k (72) 


m 


hjp — 


Kl ' 

\R> 

”55” 


■Y 1 

? > n i Z— / 


(° p ) i ni + 
R 


i=l 


m m 


4(Z - l)^ni+ 4^ni g (73) 


i=l 


i=l 


i=l 


gives the following forms for and N(k) ; 

/ m m \ 

/ E a ij a kj^ lb i £ a ij h j n j \ 

/ 0=1 0=1 \ 


N 


( 0 ) 


m 

E a kj n j 

J =1 


0 


m 

E h j“j 
0=1 


m 


i,k = 1 . . .1 (74) 


E a kj tL j n j H / RT 
0 1 


A Zt> i + 


+ S hjhjn y 

£ HjvA 


r (0) « 


m 


^ + E “J f 3 


j=l 


m 


A § + £ W; 

0=1 


i,k = 1 ... I 


I 


lA r ) = (- 1) 1 


/ a iJ n j(?f )jt a kt 0 

J,t=l 


m 


j , t=l 

m 




111 \ 

fiE^ a ij n j(lf )jt h t\ 
,t-l ^ 


0,t=l 

m 




t^t 


0 


0,t=l 


hjnjfr ^) 


CljtiH / 


r = 1,2,3 . . i,k = 1 ... I 


(75) 


(76) 
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l j a ij n j(£ r )jt f t\ 


,(*) . 


= (-i) J 


m 


j>t=l 

m 




\ <3 J t-1 


h j n j^ r } jt f t 


/ 


r = 1,2,3 . 


A considerable amount of cancellation occurs in N^ r ^ and 
(r = 1,2,3 . . .) because of two properties of the matrix T. 



n i r ik 


= 0 


(77) 


(78) 



(79) 


These properties are a direct result of Euler's theorem for homogeneous func- 
tions. The function F c , defined by equation (26), is homogeneous of degree 1 
in the np because a is homogeneous of degree zero in np and thus so is 
Z. By Euler's theorem 


m 



i=l 


^c 


(80) 


Differentiation with respect to n^ gives 


m 



i=l 


dn^dnp 


= 0 


(81) 


Equation (8l) gives immediately the second of the properties of T (eq. (79)). 
The first property (eq. (78)) also follows from this result by recalling that 
the order of differentiation is immaterial; and thus for % not equal to 
zero equation (8l) can be written in the form 

*F C 

c^iCtak 

i=l 

This is just the first property of the matrix F (eq. (78)). 
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- Equation. (78) immediately shows that the second last row of N^ r ) and the 
second last element of vanish identically. Because of equation (79) the 

contribution to N^ r ) from the term 4(Z - l) in hj vanishes as does the 
contribution to from the term 1 + 3(Z - l) - In Z in fj. 

The same expressions for N^ r ) and %( T ) (r = 0,1,2 . . .) can be used 
when the thermodynamic state is specified by temperature and pressure except 

that now the last row and last column of as well as the last element of 

^( r ) must be deleted. When entropy and pressure are used as the thermodynamic 
variables, it is only necessary to replace the last row of with 


j=l 


a kj s j n j 


S 

R 


m 


*\ 


S T + £ s j h 3 n i 

j-i 


r = 0 



s j n j(E r )jt a kt 0 


m 


s j n ^ht h t 




V k = 1,2,3 


0 > t— 1 


The last element of j £ 
A As 


(r) 


«)> t— 1 

r = 1,2,3 . 

must be replaced by 
m 


/ 


•J 


R 


+ (2 - Z - In Z)AP + 5^ sjnjfj 

J=1 


r = 0 


m 




r = 1,2,3 


j,t=l 

where the following notation was used: 


J 


(82) 


(83) 


S j = 


Klj 


bn . 

o 


+2-Z-lnZ= 


m 


T,P 


/■'"“j + 


VS 


dZ 

n i as: a j 


i=l 


Sip = T 


"5^ 


Y^(C°) m n 

\ p i 1 z + i y dz 

. - / R + Z Z~t n i da 

J n i / J 

i=l 1_x 


a/p 


(84) 


(85) 
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SUMMARY OF RESULTS 


A calculating scheme has been presented for obtaining equilibrium composi- 
tions and thermodynamic properties of plasmas containing many species. This 
scheme was made compatible with an existing calculating method for ideal gas 
mixtures by using a perturbation technique for generating a reduced set of 
working equations. 

The method permits the inclusion of both Debye-Huckel and ionization po- 
tential lowering effects in the calculation of plasma properties . These two 
effects are partially compensating. 

The equation of state was inverted in order to permit specifying the 
thermodynamic state in terms of pressure and temperature (or a function of 
temperature) rather than the usual way of assigning volume and temperature. 

The method given also permits the direct calculation of thermodynamic 
derivatives such as heat capacity and thus elimates the necessity of using the 
usually less accurate method of numerical differentiation. 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, February 17, 1965. 
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APPENDIX - SYMBOLS 


A 

A 

a ij 


a O 

°P 

E 

E^) 

0 

F 

H 

h 

fc 

I 

k 


nu 


N i 


n-? 


P 

P 

Qint 

^tr 


<li 

R 


S 


s 

T 


Helmholtz free energy- 
total mass of system 
stoichiometric coefficients 
Bohr radius 

heat capacity at constant pressure 
energy 

ground state energy level for i^ 1 isolated species 

Gibbs free energy 

enthalpy 

enthalpy per unit mass 

Planck’s constant divided by 2it 

ionization potential 

Boltzmann's constant 

mass per particle of species i 

number of particles of species i 

number of moles of species i 

absolute pressure 

partial pressure 

internal partition function 

translational partition function 

charge on i^ species 

universal gas constant 

entropy 

entropy per unit mass 
absolute temperature 



U internal energy 

V volume 

W mass 

Z Coulomb compressibility 

a dimensionless ionization parameter, eq. (27) 

3j_ proportionality constant 

T perturbation matrix, eq. (43) 

6 . . Kronecker delta 

**• J 

isolated energy level relative to ground state 
k reciprocal of Debye length 

« Lagrange multiplier 

Subscripts : 

c Coulomb forces 

I ideal gas 

0 ground state 

Superscripts : 
o standard state 
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ERRATA. 


NASA Technical Note D-2806 

EQUILIBRIUM COMPUTATIONS FOR MULTICOMPONENT PLASMAS 

By Frank J. Zeleznik and Sanford Gordon 
May 1965 

Page 15: Line 8 should read 

M£l = (I + A) -1 , where A = r except for equations (57) and 
(58). The need for numerically calculating would be 

Page 15: All T's below line 8 should be changed to A's. 

Page 16, equation (67): The P's should be changed to A*s. 
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